SUBROUTINE sample(element,s,wt)
!
! This subroutine returns the local coordinates and weighting coefficients
! of the integrating points.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(OUT)::s(:,:)
 REAL(iwp),INTENT(OUT),OPTIONAL::wt(:)
 CHARACTER(*),INTENT(IN)::element
 INTEGER::nip
 REAL(iwp)::root3,r15,w(3),v(9),b,c
 root3=1.0_iwp/SQRT(3.0_iwp)
 r15=0.2_iwp*SQRT(15.0_iwp)
 nip=UBOUND(s,1)
 w=(/5.0_iwp/9.0_iwp,8.0_iwp/9.0_iwp,5.0_iwp/9.0_iwp/)
 v=(/5.0_iwp/9.0_iwp*w,8.0_iwp/9.0_iwp*w,5.0_iwp/9.0_iwp*w/)
 SELECT CASE(element)
 CASE('line')
   SELECT CASE(nip)
   CASE(1)
     s(1,1)=0.0_iwp
     wt(1) =2.0_iwp
   CASE(2)
     s(1,1)=-0.577350269189626_iwp
     s(2,1)= 0.577350269189626_iwp
     wt(1) = 1.000000000000000_iwp
     wt(2) = 1.000000000000000_iwp
   CASE(3)
     s(1,1)=-0.774596669241484_iwp
     s(2,1)= 0.000000000000000_iwp
     s(3,1)= 0.774596669241484_iwp
     wt(1) = 0.555555555555556_iwp
     wt(2) = 0.888888888888889_iwp
     wt(3) = 0.555555555555556_iwp
   CASE(4)
     s(1,1)=-0.861136311594053_iwp
     s(2,1)=-0.339981043584856_iwp
     s(3,1)= 0.339981043584856_iwp
     s(4,1)= 0.861136311594053_iwp
     wt(1) = 0.347854845137454_iwp
     wt(2) = 0.652145154862546_iwp
     wt(3) = 0.652145154862546_iwp
     wt(4) = 0.347854845137454_iwp
   CASE(5)
     s(1,1)=-0.906179845938664_iwp
     s(2,1)=-0.538469310105683_iwp
     s(3,1)= 0.000000000000000_iwp
     s(4,1)= 0.538469310105683_iwp
     s(5,1)= 0.906179845938664_iwp
     wt(1) = 0.236926885056189_iwp
     wt(2) = 0.478628670499366_iwp
     wt(3) = 0.568888888888889_iwp
     wt(4) = 0.478628670499366_iwp
     wt(5) = 0.236926885056189_iwp
   CASE(6)
     s(1,1)=-0.932469514203152_iwp
     s(2,1)=-0.661209386466265_iwp
     s(3,1)=-0.238619186083197_iwp
     s(4,1)= 0.238619186083197_iwp
     s(5,1)= 0.661209386466265_iwp
     s(6,1)= 0.932469514203152_iwp
     wt(1) = 0.171324492379170_iwp
     wt(2) = 0.360761573048139_iwp
     wt(3) = 0.467913934572691_iwp
     wt(4) = 0.467913934572691_iwp
     wt(5) = 0.360761573048139_iwp
     wt(6) = 0.171324492379170_iwp
   CASE(7)
     s(1,1)=-0.9491079123427585245261897_iwp
     s(2,1)=-0.7415311855993944398638648_iwp
     s(3,1)=-0.4058451513773971669066064_iwp
     s(4,1)= 0.000000000000000_iwp
     s(5,1)= 0.4058451513773971669066064_iwp
     s(6,1)= 0.7415311855993944398638648_iwp
     s(7,1)= 0.9491079123427585245261897_iwp
     wt(1) = 0.1294849661688696932706114_iwp
     wt(2) = 0.2797053914892766679014678_iwp
     wt(3) = 0.3818300505051189449503698_iwp
     wt(4) = 0.4179591836734693877551020_iwp
     wt(5) = 0.3818300505051189449503698_iwp
     wt(6) = 0.2797053914892766679014678_iwp
     wt(7) = 0.1294849661688696932706114_iwp
   CASE(8)
     s(1,1)=-0.9602898564975362316835609_iwp
     s(2,1)=-0.7966664774136267395915539_iwp
     s(3,1)=-0.5255324099163289858177390_iwp
     s(4,1)=-0.1834346424956498049394761_iwp
     s(5,1)= 0.1834346424956498049394761_iwp
     s(6,1)= 0.5255324099163289858177390_iwp
     s(7,1)= 0.7966664774136267395915539_iwp
     s(8,1)= 0.9602898564975362316835609_iwp
     wt(1) = 0.1012285362903762591525314_iwp
     wt(2) = 0.2223810344533744705443560_iwp
     wt(3) = 0.3137066458778872873379622_iwp
     wt(4) = 0.3626837833783619829651504_iwp
     wt(5) = 0.3626837833783619829651504_iwp
     wt(6) = 0.3137066458778872873379622_iwp
     wt(7) = 0.2223810344533744705443560_iwp
     wt(8) = 0.1012285362903762591525314_iwp
   CASE(9)
     s(1,1)=-0.9681602395076260898355762_iwp
     s(2,1)=-0.8360311073266357942994298_iwp    
     s(3,1)=-0.6133714327005903973087020_iwp
     s(4,1)=-0.3242534234038089290385380_iwp    
     s(5,1)= 0.000000000000000_iwp                            
     s(6,1)= 0.3242534234038089290385380_iwp                            
     s(7,1)= 0.6133714327005903973087020_iwp                            
     s(8,1)= 0.8360311073266357942994298_iwp                            
     s(9,1)= 0.9681602395076260898355762_iwp                            
     wt(1) = 0.0812743883615744119718922_iwp                            
     wt(2) = 0.1806481606948574040584720_iwp                            
     wt(3) = 0.2606106964029354623187429_iwp                            
     wt(4) = 0.3123470770400028400686304_iwp                            
     wt(5) = 0.3302393550012597631645251_iwp                            
     wt(6) = 0.3123470770400028400686304_iwp                            
     wt(7) = 0.2606106964029354623187429_iwp                            
     wt(8) = 0.1806481606948574040584720_iwp                            
     wt(9) = 0.0812743883615744119718922_iwp                            
   CASE(10)
     s(1,1)=-0.9739065285171717200779640_iwp            
     s(2,1)=-0.8650633666889845107320967_iwp 
     s(3,1)=-0.6794095682990244062343274_iwp 
     s(4,1)=-0.4333953941292471907992659_iwp 
     s(5,1)=-0.1488743389816312108848260_iwp 
     s(6,1)= 0.1488743389816312108848260_iwp 
     s(7,1)= 0.4333953941292471907992659_iwp 
     s(8,1)= 0.6794095682990244062343274_iwp 
     s(9,1)= 0.8650633666889845107320967_iwp 
    s(10,1)= 0.9739065285171717200779640_iwp 
     wt(1) = 0.0666713443086881375935688_iwp                     
     wt(2) = 0.1494513491505805931457763_iwp                     
     wt(3) = 0.2190863625159820439955349_iwp                     
     wt(4) = 0.2692667193099963550912269_iwp                     
     wt(5) = 0.2955242247147528701738930_iwp                     
     wt(6) = 0.2955242247147528701738930_iwp                      
     wt(7) = 0.2692667193099963550912269_iwp                     
     wt(8) = 0.2190863625159820439955349_iwp                     
     wt(9) = 0.1494513491505805931457763_iwp                     
    wt(10) = 0.0666713443086881375935688_iwp                     
   CASE DEFAULT                              
     WRITE(*,*)"Wrong number of integrating points for a line"
   END SELECT
 CASE('triangle')
 SELECT CASE(nip)
   CASE(1)
     s(1,1)= 0.333333333333333_iwp
     s(1,2)= 0.333333333333333_iwp
     wt(1) = 0.500000000000000_iwp
   CASE(3)
     s(1,1)= 0.500000000000000_iwp
     s(1,2)= 0.500000000000000_iwp
     s(2,1)= 0.500000000000000_iwp
     s(2,2)= 0.000000000000000_iwp
     s(3,1)= 0.000000000000000_iwp
     s(3,2)= 0.500000000000000_iwp
     wt(1:3)=0.333333333333333_iwp
     wt=0.5_iwp*wt
   CASE(4)
     s(1,1)= 0.6_iwp
     s(1,2)= 0.2_iwp
     s(2,1)= 0.2_iwp
     s(2,2)= 0.6_iwp
     s(3,1)= 0.2_iwp
     s(3,2)= 0.2_iwp
     s(4,1)= 0.333333333333333_iwp
     s(4,2)= 0.333333333333333_iwp
     wt(1:3)= 0.520833333333333_iwp
     wt(4)=  -0.5625_iwp
     wt=0.5_iwp*wt
   CASE(6)
     s(1,1)= 0.816847572980459_iwp
     s(1,2)= 0.091576213509771_iwp
     s(2,1)= 0.091576213509771_iwp
     s(2,2)= 0.816847572980459_iwp
     s(3,1)= 0.091576213509771_iwp
     s(3,2)= 0.091576213509771_iwp
     s(4,1)= 0.108103018168070_iwp
     s(4,2)= 0.445948490915965_iwp
     s(5,1)= 0.445948490915965_iwp
     s(5,2)= 0.108103018168070_iwp
     s(6,1)= 0.445948490915965_iwp
     s(6,2)= 0.445948490915965_iwp
     wt(1:3)=0.109951743655322_iwp
     wt(4:6)=0.223381589678011_iwp
     wt=0.5_iwp*wt
   CASE(7)
     s(1,1)= 0.333333333333333_iwp
     s(1,2)= 0.333333333333333_iwp
     s(2,1)= 0.797426985353087_iwp
     s(2,2)= 0.101286507323456_iwp
     s(3,1)= 0.101286507323456_iwp
     s(3,2)= 0.797426985353087_iwp
     s(4,1)= 0.101286507323456_iwp
     s(4,2)= 0.101286507323456_iwp
     s(5,1)= 0.470142064105115_iwp
     s(5,2)= 0.059715871789770_iwp
     s(6,1)= 0.059715871789770_iwp
     s(6,2)= 0.470142064105115_iwp
     s(7,1)= 0.470142064105115_iwp
     s(7,2)= 0.470142064105115_iwp
     wt(1) = 0.225000000000000_iwp
     wt(2:4)=0.125939180544827_iwp
     wt(5:7)=0.132394152788506_iwp
     wt=0.5_iwp*wt
   CASE(12)
     s(1,1)= 0.873821971016996_iwp
     s(1,2)= 0.063089014491502_iwp
     s(2,1)= 0.063089014491502_iwp
     s(2,2)= 0.873821971016996_iwp
     s(3,1)= 0.063089014491502_iwp
     s(3,2)= 0.063089014491502_iwp
     s(4,1)= 0.501426509658179_iwp
     s(4,2)= 0.249286745170910_iwp
     s(5,1)= 0.249286745170910_iwp
     s(5,2)= 0.501426509658179_iwp
     s(6,1)= 0.249286745170910_iwp
     s(6,2)= 0.249286745170910_iwp
     s(7,1) =0.053145049844817_iwp
     s(7,2) =0.310352451033784_iwp
     s(8,1) =0.310352451033784_iwp
     s(8,2) =0.053145049844817_iwp
     s(9,1) =0.053145049844817_iwp
     s(9,2) =0.636502499121398_iwp
     s(10,1)=0.310352451033784_iwp
     s(10,2)=0.636502499121398_iwp
     s(11,1)=0.636502499121398_iwp
     s(11,2)=0.053145049844817_iwp
     s(12,1)=0.636502499121398_iwp
     s(12,2)=0.310352451033784_iwp
     wt(1:3)=0.050844906370207_iwp
     wt(4:6)=0.116786275726379_iwp
     wt(7:12)=0.082851075618374_iwp
     wt=0.5_iwp*wt
   CASE(16)
     s(1,1)=0.333333333333333_iwp
     s(1,2)=0.333333333333333_iwp
     s(2,1)=0.658861384496478_iwp
     s(2,2)=0.170569307751761_iwp
     s(3,1)=0.170569307751761_iwp
     s(3,2)=0.658861384496478_iwp
     s(4,1)=0.170569307751761_iwp
     s(4,2)=0.170569307751761_iwp
     s(5,1)=0.898905543365938_iwp
     s(5,2)=0.050547228317031_iwp
     s(6,1)=0.050547228317031_iwp
     s(6,2)=0.898905543365938_iwp
     s(7,1)=0.050547228317031_iwp
     s(7,2)=0.050547228317031_iwp
     s(8,1)=0.081414823414554_iwp
     s(8,2)=0.459292588292723_iwp
     s(9,1)=0.459292588292723_iwp
     s(9,2)=0.081414823414554_iwp
     s(10,1)=0.459292588292723_iwp
     s(10,2)=0.459292588292723_iwp
     s(11,1)=0.008394777409958_iwp
     s(11,2)=0.263112829634638_iwp
     s(12,1)=0.008394777409958_iwp
     s(12,2)=0.728492392955404_iwp
     s(13,1)=0.263112829634638_iwp
     s(13,2)=0.008394777409958_iwp
     s(14,1)=0.263112829634638_iwp
     s(14,2)=0.728492392955404_iwp
     s(15,1)=0.728492392955404_iwp
     s(15,2)=0.008394777409958_iwp
     s(16,1)=0.728492392955404_iwp
     s(16,2)=0.263112829634638_iwp
     wt(1)=0.144315607677787_iwp
     wt(2:4)=0.103217370534718_iwp
     wt(5:7)=0.032458497623198_iwp
     wt(8:10)=0.095091634267284_iwp
     wt(11:16)=0.027230314174435_iwp
     wt=0.5_iwp*wt
   CASE DEFAULT
     WRITE(*,*)"wrong number of integrating points for a triangle"
   END SELECT
 CASE('quadrilateral')
   SELECT CASE(nip)
   CASE(1)
     s(1,1)=0.0_iwp
     s(1,2)=0.0_iwp
     wt(1)=4.0_iwp
   CASE(4)
     s(1,1)=-root3
     s(1,2)= root3
     s(2,1)= root3
     s(2,2)= root3
     s(3,1)=-root3
     s(3,2)=-root3
     s(4,1)= root3
     s(4,2)=-root3
     wt=1.0_iwp
   CASE(9)
     s(1:7:3,1)=-r15
     s(2:8:3,1)=0.0_iwp
     s(3:9:3,1)=r15
     s(1:3,2)  =r15
     s(4:6,2)  =0.0_iwp
     s(7:9,2)  =-r15
     wt= v
   CASE(16)
     s(1:13:4,1)=-0.861136311594053_iwp
     s(2:14:4,1)=-0.339981043584856_iwp
     s(3:15:4,1)= 0.339981043584856_iwp
     s(4:16:4,1)= 0.861136311594053_iwp
     s(1:4,2)   = 0.861136311594053_iwp
     s(5:8,2)   = 0.339981043584856_iwp
     s(9:12,2)  =-0.339981043584856_iwp
     s(13:16,2) =-0.861136311594053_iwp
     wt(1)      = 0.121002993285602_iwp
     wt(4)      = wt(1)
     wt(13)     = wt(1)
     wt(16)     = wt(1)
     wt(2)      = 0.226851851851852_iwp
     wt(3)      = wt(2)
     wt(5)      = wt(2)
     wt(8)      = wt(2)
     wt(9)      = wt(2)
     wt(12)     = wt(2)
     wt(14)     = wt(2)
     wt(15)     = wt(2)
     wt(6)      = 0.425293303010694_iwp
     wt(7)      = wt(6)
     wt(10)     = wt(6)
     wt(11)     = wt(6)
   CASE(25)
     s(1:21:5,1)= 0.906179845938664_iwp
     s(2:22:5,1)= 0.538469310105683_iwp
     s(3:23:5,1)= 0.0_iwp
     s(4:24:5,1)=-0.538469310105683_iwp
     s(5:25:5,1)=-0.906179845938664_iwp
     s( 1: 5,2) = 0.906179845938664_iwp
     s( 6:10,2) = 0.538469310105683_iwp
     s(11:15,2) = 0.0_iwp
     s(16:20,2) =-0.538469310105683_iwp
     s(21:25,2) =-0.906179845938664_iwp
     wt(1) =0.056134348862429_iwp
     wt(2) =0.113400000000000_iwp
     wt(3) =0.134785072387521_iwp
     wt(4) =0.113400000000000_iwp
     wt(5) =0.056134348862429_iwp
     wt(6) =0.113400000000000_iwp
     wt(7) =0.229085404223991_iwp
     wt(8) =0.272286532550750_iwp
     wt(9) =0.229085404223991_iwp
     wt(10)=0.113400000000000_iwp
     wt(11)=0.134785072387521_iwp
     wt(12)=0.272286532550750_iwp
     wt(13)=0.323634567901235_iwp
     wt(14)=0.272286532550750_iwp
     wt(15)=0.134785072387521_iwp
     wt(16)=0.113400000000000_iwp
     wt(17)=0.229085404223991_iwp
     wt(18)=0.272286532550750_iwp
     wt(19)=0.229085404223991_iwp
     wt(20)=0.113400000000000_iwp
     wt(21)=0.056134348862429_iwp
     wt(22)=0.113400000000000_iwp
     wt(23)=0.134785072387521_iwp
     wt(24)=0.113400000000000_iwp
     wt(25)=0.056134348862429_iwp
   CASE DEFAULT
     WRITE(*,*)"wrong number of integrating points for a quadrilateral"
   END SELECT
 CASE('tetrahedron')
!                       for tetrahedra weights multiplied by 1/6
   SELECT CASE(nip)
   CASE(1)
     s(1,1)=0.25_iwp
     s(1,2)=0.25_iwp
     s(1,3)=0.25_iwp
     wt(1)=1.0_iwp/6.0_iwp
   CASE(4)
     s(1,1)=0.58541020_iwp
     s(1,2)=0.13819660_iwp
     s(1,3)=s(1,2)
     s(2,2)=s(1,1)
     s(2,3)=s(1,2)
     s(2,1)=s(1,2)
     s(3,3)=s(1,1)
     s(3,1)=s(1,2)
     s(3,2)=s(1,2)
     s(4,1)=s(1,2)
     s(4,2)=s(1,2)
     s(4,3)=s(1,2)
     wt(1:4)=0.25_iwp/6.0_iwp
   CASE(5)
     s(1,1)=0.25_iwp
     s(1,2)=0.25_iwp
     s(1,3)=0.25_iwp
     s(2,1)=0.5_iwp
     s(2,2)=1.0_iwp/6.0_iwp
     s(2,3)=s(2,2)
     s(3,2)=0.5_iwp
     s(3,3)=1.0_iwp/6.0_iwp
     s(3,1)=s(3,3)
     s(4,3)=0.5_iwp
     s(4,1)=1.0_iwp/6.0_iwp
     s(4,2)=s(4,1)
     s(5,1)=1.0_iwp/6.0_iwp
     s(5,2)=s(5,1)
     s(5,3)=s(5,1)
     wt(1)=-0.8_iwp
     wt(2)=9.0_iwp/20.0_iwp
     wt(3:5)=wt(2)
     wt=wt/6.0_iwp
   CASE DEFAULT
     WRITE(*,*)"wrong number of integrating points for a tetrahedron"
   END SELECT
 CASE('hexahedron')
   SELECT CASE(nip)
   CASE(1)
     s(1,1:3)=0.0_iwp
     wt(1)=8.0_iwp
   CASE(8)
     s(1,1)= root3
     s(1,2)= root3
     s(1,3)= root3
     s(2,1)= root3
     s(2,2)= root3
     s(2,3)=-root3
     s(3,1)= root3
     s(3,2)=-root3
     s(3,3)= root3
     s(4,1)= root3
     s(4,2)=-root3
     s(4,3)=-root3
     s(5,1)=-root3
     s(5,2)= root3
     s(5,3)= root3
     s(6,1)=-root3
     s(6,2)=-root3
     s(6,3)= root3
     s(7,1)=-root3
     s(7,2)= root3
     s(7,3)=-root3
     s(8,1)=-root3
     s(8,2)=-root3
     s(8,3)=-root3
     wt=1.0_iwp
   CASE(14)
     b=0.795822426_iwp
     c=0.758786911_iwp
     wt(1:6)=0.886426593_iwp
     wt(7:14)=0.335180055_iwp
     s(1,1)=-b
     s(2,1)=b
     s(3,2)=-b
     s(4,2)=b
     s(5,3)=-b
     s(6,3)=b
     s(7:,:)=c
     s(7,1)=-c
     s(7,2)=-c
     s(7,3)=-c
     s(8,2)=-c
     s(8,3)=-c
     s(9,1)=-c
     s(9,3)=-c
     s(10,3)=-c
     s(11,1)=-c
     s(11,2)=-c
     s(12,2)=-c
     s(13,1)=-c
   CASE(15)
     b=1.0_iwp
     c       =0.674199862_iwp
     wt(1)   =1.564444444_iwp
     wt(2:7) =0.355555556_iwp
     wt(8:15)=0.537777778_iwp
     s(2,1)=-b
     s(3,1)=b
     s(4,2)=-b
     s(5,2)=b
     s(6,3)=-b
     s(7,3)=b
     s(8:,:)=c
     s(8,1)=-c
     s(8,2)=-c
     s(8,3)=-c
     s(9,2)=-c
     s(9,3)=-c
     s(10,1)=-c
     s(10,3)=-c
     s(11,3)=-c
     s(12,1)=-c
     s(12,2)=-c
     s(13,2)=-c
     s(14,1)=-c
   CASE(27)
     wt=(/5.0_iwp/9.0_iwp*v,8.0_iwp/9.0_iwp*v,5.0_iwp/9.0_iwp*v/)
     s(1:7:3,1)=-r15
     s(2:8:3,1)=0.0_iwp
     s(3:9:3,1)=r15
     s(1:3,3)=r15
     s(4:6,3)=0.0_iwp
     s(7:9,3)=-r15
     s(1:9,2)=-r15
     s(10:16:3,1)=-r15
     s(11:17:3,1)=0.0_iwp
     s(12:18:3,1)=r15
     s(10:12,3)=r15
     s(13:15,3)=0.0_iwp
     s(16:18,3)=-r15
     s(10:18,2)=0.0_iwp
     s(19:25:3,1)=-r15
     s(20:26:3,1)=0.0_iwp
     s(21:27:3,1)=r15
     s(19:21,3)= r15
     s(22:24,3)=0.0_iwp
     s(25:27,3)=-r15
     s(19:27,2)= r15
   CASE DEFAULT
     WRITE(*,*)"wrong number of integrating points for a hexahedron"
   END SELECT
 CASE DEFAULT
   WRITE(*,*)"not a valid element type"
 END SELECT
RETURN
END SUBROUTINE sample
